cc     Program by Pablo I. Hurtado, based on code by Walter Kob for a Lennard-Jones binary mixture.
c       polydisperse softspheres



      implicit double precision (a-h,o-z)

c     Parameters: 

c     ipartdim1: number of particles dimension (to define vectors)

c     ipartdim: number of particle coordinates (3*ipartdim1)

c     ipoldim: number of polymers dimension

c     ipoldim2: number of polymer endpoints dimension

c     idimen: second dimension for the polymer-induced neighbor list

      parameter (ipartdim1=10000,idimen=600,idimen2=150)

      parameter (ipartdim=3*ipartdim1,ipoldim=20*ipartdim1)

      parameter(ipoldim2=2*ipoldim)

      parameter(icorrdim=400,icorrdim2=10000)



      dimension kseed(24),neighboro(idimen2*ipartdim1),

     &numpol(ipartdim1),neighpol(ipartdim1,idimen),

     &kpol(ipartdim1,idimen),kendpoint(ipoldim2),indexpol(ipoldim),

     &nvecp(ipartdim1),kvecpolymer(ipartdim1,idimen)

      dimension pos(ipartdim),vel(ipartdim),acc(ipartdim),

     &runpos(ipartdim)

      dimension  histox(-100:100),histoy(-100:100),histoz(-100:100),

     &histor(0:100),histovx(-100:100),histovy(-100:100),

     &histovz(-100:100),histov(0:100)

      dimension icortim(icorrdim),icortimwt(icorrdim2)

      dimension itimewrite(icorrdim*icorrdim2)



      character*200 descrip,descrip0

      character*300 filein0,filein1,filein2,fileout1,fileout2,

     &fileout4,filein3,fichier

      character*300 dump_file,fileout3,filepot

      character*6 cnpart

      character*9 cfpol,campli,cL,cvfrac,ctemp,cstep,ctpol,cftrial

c      character*132 snap_file,snap_file1,string1

c --------------------------------------------------------------

       real*8 dia(ipartdim1),dia3

      common/kob1/ xlbox,xlboxd2,npart,npol,ntot

      common/kob3/ hstep,ekin,epot,epotpol,fracloop

      common/kob4/ rcutoff, rcutoff2

      common/kob7/ rskin2

      common/kob8/ vel

      common/kob9/ tdtp24,icarry,i24,j24,kseed,itwop24

      common/kob10/ neighboro                    ! neighboro is the old neighbor list; neighbor is the new one

      common/kob11/ virial

      common/pablo1/ npart3

      common/pablo2/ neighpol,numpol

      common/pablo3/ cte1,cte2,cte3,cte4,polmaxL,polmaxL2,Rcore,epsilon

      common/pablo4/ ctep1,ctep2,ctep3,Rcore2,Rcore14

      common/pablo5/ distmax2,kcountmax

      common/pablo6/ nvecp,kvecpolymer

      common/pablo7/ aparam,bparam,cparam,fLL,fLL2

      common/pablo8/ hkk,fmax1,fmax2,distm1,distm2,itime

      common/measure1/ histox,histoy,histoz,histor

      common/measure2/ histovx,histovy,histovz,histov

      common/mpara/ vmax,dhv,dhx,xpmed,npoints

      common/pc1/dia


c      open(066,file='input.dat',status='old')

c      read(066,fmt='(a)') filein0

c      close(066)



c      open(068,file=filein0,status='old')

c      read(068,fmt='(a)') filein1

c      read(068,fmt='(a)') fileout1

c      read(068,fmt='(a)') fileout2

c      read(068,fmt='(a)') fileout3

c      read(068,fmt='(a)') fileout4

c      read(068,fmt='(a)') filein2

c      read(068,fmt='(a)') filepot

c      read(068,fmt='(a)') filein3

c      close(068)



      kcountmax=idimen2*ipartdim1



c

c read in the start parameter file

c



      open(061,file='xstartparameter-file.dat',status='old')

      read(061,*) npart                          ! Number of particles

      read(061,*) npol                            ! Number of polymers

c     data for hard-core and polymer-induced potentials

      read(061,*) Rcore                        ! Length scale for the hard-core interaction (~ particle DIAMETER)

      read(061,*) rcutoff                        ! Cutoff of "free" potential (for disconnected particles)

      read(061,*) ampli                         ! Amplitude of the polymer-mediated potential

      read(061,*) polmaxL                     ! Polymer maximum length

      read(061,*) polR0                         ! Polymer radius of giration

      read(061,*) epsilon                      ! Free energy cost of a looped polymer

c     system data

      read(061,*) volfrac                       ! Particle volume fraction

      read(061,*) rskin                          ! Skin for the Verlet list

      read(061,*) temp                          ! Temperature

      read(061,*) icooltim                       ! Time to thermalize the system

      read(061,*) icoolstep                    ! Time interval between two succesive velocity randomizations

      read(061,*) etot                            ! Total energy PER PARTICLE

      read(061,*) iolseed                       ! Old seed for random number gnerator

      read(061,*) hstep                         ! Time step of Verlet MD scheme

      read(061,*) ntotstep                     ! Total number of time steps

      read(061,*) nmedida                    ! Time interval between succesive measurements

      read(061,*) nprintstep                  ! Time interval between succesive printings to file: SHOULD BE MULTIPLE OF nmedida

      read(061,*) ndumpstep                ! Time step between succesive info dumpings to file

      read(061,*) itneighli                      ! Time step between actualizations of Verlet list

      read(061,*) timepol                       ! Average time between different polymer network rearrangements

      read(061,*) fractrial                      ! Fraction of Monte Carlo move trials for polymers at each Verlet list update

      read(061,*) fL                               ! Fraction of maximum polymer extension we use to regularize the FENE

c                                                         ! force law in order to avoid spurious super-extensions.

      read(061,*) ncorrel                      ! Number of log-points for the correlation functions, mean-square displacement, etc

      read(061,*) ncorrelwt                   ! Number of linear waiting times for correlation functions

      close(061)



      write(cnpart,fmt='(i6.6)')npart

c      write(cnpol,fmt='(i6.6)')npol

      write(cfpol,fmt='(e9.4)')2.*npol/npart

      write(campli,fmt='(e9.4)')ampli

      write(cL,fmt='(e9.4)')polmaxL

      write(cvfrac,fmt='(e9.4)')volfrac

      write(ctemp,fmt='(e9.4)')temp

      write(cstep,fmt='(e9.4)')hstep

      write(ctpol,fmt='(e9.4)')timepol

      write(cftrial,fmt='(e9.4)')fractrial




      descrip0='-N_'//cnpart//'-vol_'//cvfrac//'-rpol_'//cfpol//'-T_'/

     &/ctemp//'-A_'//campli//'-Lmax_'//cL//'-h_'//cstep//'-tau_'/

     &/ctpol//'-trial_'//cftrial//'.'



      descrip='-N_'//cnpart//'-vol_'//cvfrac//'-rpol_'//cfpol//'-T_'/

     &/ctemp//'-A_'//campli//'-Lmax_'//cL//'-h_'//cstep//'-tau_'/

     &/ctpol//'-trial_'//cftrial//'.sal'



      fileout1='log-file'//descrip

      fileout2='dump-file'//descrip

      fileout3='final-config'//descrip

      fileout4='snapshot-file'//descrip

      filepot='potential-force'//descrip

      

      filein0='config'//descrip0

      filein2='start-config.dat'

      filein3='yhistograms.dat'

c ---------------------------------------------------------------------------------

	pi=0.0d0
	pi=4.0d0*datan(1.0d0)


        d=0.150d0

	a=0.0d0
	b=0.0d0
	dia3=0.0d0
	dia2=0.0d0
	dia=0.0d0

        b=1+d*sqrt(3.0)
        a=1-d*sqrt(3.0)

	open(unit=45,file='dia.dat',status='unknown')

        do i=1,npart
c        dia(i)=a+(b-a)*rand()
        read(45,*)dia(i)
        dia3=dia3+dia(i)**3.0
	dia1=dia1+dia(i)
	dia2=dia2+dia(i)**2.0d0
        enddo

	close(45)

	dia3=dia3/dfloat(npart)
	dia1=dia1/dfloat(npart)
	dia2=dia2/dfloat(npart)

	print*,dia3,dia1,dsqrt(dia2/(dia1*dia1)-1.0d0)


c -------------------------------------------------------------------------------



     
c        xlbox=Rcore*(dacos(-1.d0)*npart/(6.d0*volfrac))**(1.d0/3.d0) 

c       print*,xlbox

	xlbox=0.0d0

         xlbox=(dia3*dacos(-1.d0)*npart/(6.d0*volfrac))**(1.d0/3.d0)          ! System linear size

         print*,xlbox,(pi/6.0d0)*(dfloat(npart)/(xlbox)**3.0d0)*dia3

	 gamm=0.0d0
	 gamm=(dfloat(npart)/(xlbox)**3.0d0)*dia3
	 gamm=gamm*(temp**(-3.0d0/14.0d0))

	print*,d,gamm


 
c        xlbox=Rcore*(dacos(-1.d0)*npart/(6.d0*volfrac))**(1.d0/3.d0) 



      ntrial=idint(fractrial*npol)                           ! Number of polymers we try to move at every Verlet list update

      xnu0=hstep*itneighli/(fractrial*timepol)        ! Intrinsic acceptance frecuency for Metropolis rate

      xcte=xnu0*timepol

      if(xnu0.gt.1.d0) then

         write(*,*)'xnu0>1   ERROR'

         stop

      endif



      open(062,file=filein3,status='old')

      read(062,*) npoints                               ! Number of points in velocity and spatial histograms: SHOULD BE EVEN

      read(062,*) vmax                                  ! Maximum velocity to be considered in the velocity histograms

      close(62)



      dhv=vmax/npoints

      dhx=(Rcore+polmaxL)/npoints

      xpmed=0

      xmed=0

      do i=0,100

         histox(i)=0

         histox(-i)=0

         histoy(i)=0

         histoy(-i)=0

         histoz(i)=0

         histoz(-i)=0



         histovx(i)=0

         histovx(-i)=0

         histovy(i)=0

         histovy(-i)=0

         histovz(i)=0

         histovz(-i)=0



         histor(i)=0

         histov(i)=0

      enddo

      

      open(071,file=fileout1,status='unknown',form='formatted')

      write(071,*)'#  ',npart,'    npart'

      write(071,*)'#  ',npol,'    npol'

      write(071,*)'#  ',Rcore,'    Rcore'

      write(071,*)'#  ',rcutoff,'    rcutoff'

      write(071,*)'#  ',ampli,'    ampli'

      write(071,*)'#  ',polmaxL,'   polmaxL'

      write(071,*)'#  ',polR0,'    polR0'

      write(071,*)'#  ',epsilon,'    epsilon'

      write(071,*)'#  ',volfrac,'    volfrac'

      write(071,*)'#  ',xlbox,'    xlbox'

      write(071,*)'#  ',rskin,'    rskin'

      write(071,*)'#  ',temp,'    temp'

      write(071,*)'#  ',icooltim,'    icooltim'

      write(071,*)'#  ',icoolstep,'    icoolstep'

      write(071,*)'#  ',etot,'    etot'

      write(071,*)'#  ',iolseed,'    iolseed'

      write(071,*)'#  ',hstep,'    hstep'

      write(071,*)'#  ',ntotstep,'    ntotstep'

      write(071,*)'#  ',nmedida,'    nmedida'

      write(071,*)'#  ',nprintstep,'    nprintstep'

      write(071,*)'#  ',ndumpstep,'    ndumpstep'

      write(071,*)'#  ',itneighli,'    itneighli'

      write(071,*)'#  ',timepol,'    timepol'

      write(071,*)'#  ',fractrial,'    fractrial'

      write(071,*)'#  ',ntrial,'    ntrial'

      write(071,*)'#  ',xnu0,'    xnu0'



c     Compute the log-times when the correlation functions are evaluated. The time

c     interval over which correlation functions are measured corresponds to

c     ndumpstep molecular dynamics steps.

      do k=1,ncorrel

c         icortim(k)=nint(dfloat(ndumpstep)**(k/dfloat(ncorrel+1)))

         icortim(k)=nint(dfloat(ndumpstep)**(k/dfloat(ncorrel)))

         write(*,*)icortim(k)

      enddo

      icortim(1)=0



c      pause



c     Since the vector icortim(k) is the integer part of a real number, there exists the 

c     possibility that two consecutive vectors cortim(k) and cortim(k+1) have the same 

c     integer value. In order to correct this problem, if present, we perform the following loop.

      k=1

      j=2

 20   if( j .gt. ncorrel ) goto 30

      if( icortim(j) .ne. icortim(k) ) then

         k=k+1

         icortim(k)=icortim(j)

         j=j+1

      else

         j=j+1

      endif

      goto 20



 30   if( ncorrel .eq. 0 .or. ndumpstep .le. 1 ) then

         ncorrelnew=0

      else

         ncorrelnew=k

      endif



      write(071,*) '# log-times'

      write(071,*) '# ncorrelnew=',ncorrelnew

      do k=1,ncorrelnew

         write(071,*) '#    ',icortim(k)

c         write(33,*)icortim(k)

      enddo

c      call flushkob(33)


c      if( ncorrelnew+ncorrelmax+1 .gt. cor1dim+cor2dim ) then

c        print*,'error5', ncorrel,ncorrelnew,ncorrelmax,cor1dim,cor2dim

c        write(071,*) 'error5', ncorrel,ncorrelnew,ncorrelmax,cor1dim,

c     &                cor2dim

c        stop

c      endif



c     compute the times of the linear waiting times

      j=ncorrelwt

      ncorrelwtnew=ncorrelwt

      do k=1,j

         icortimwt(k)=(k-1)*nint(ntotstep/dfloat(j))

      enddo

      write(071,*) '# linear waiting times'

      write(071,*) '# ncorrelwtnew=',ncorrelwtnew

      ii=0

      do k=1,ncorrelwtnew

         write(071,*)  '#    ',icortimwt(k)

c         write(34,*)icortimwt(k)

         do l=1,ncorrelnew

            ii=ii+1

            itimewrite(ii)=icortim(l)+icortimwt(k)

         enddo

      enddo

c      call flushkob(34)



c     We now order the writing times list itimewrite(ii), and initialize

c     the writing counter nwrite to 1. The number iwr is the total number

c     of writing times.

      iwr=ncorrelnew*ncorrelwtnew                   ! Total number of times for which we write the configuration

c      do ii=1,iwr

c         write(35,*)itimewrite(ii)

c      enddo

c      call flushkob(35)

      call sort(iwr,itimewrite)

      nwrite=1

      do ii=1,iwr

         write(36,*)itimewrite(ii)

      enddo

      call flushkob(36)



c      pause



c

c set up some constants

c



      hkk=hstep

      xaccept=0

      xmov=0

      fmax1=0.d0

      fmax2=0.d0

      distm1=100.d0

      distm2=0.d0



c for the rng



      itwop24=16777216

      xtwop24=1.d0/dfloat(itwop24)

      i24=24

      j24=10

      icarry=0



      tdtp24=2.0/dfloat(itwop24+1)



c for the simulation



      ntot=npart+npol

      npart3=3*npart

      density=npart/xlbox**3

      rcutoff2=rcutoff*rcutoff

      rcutoff4=rcutoff2*rcutoff2

      npol2=2*npol



      beta=1.0/temp



      rskin2=rcutoff2*rskin**2

      xlboxd2=xlbox/2.d0

      hsqd2=hstep**2/2.d0

      hstepd2=hstep/2.d0



c      Rcore2=Rcore*Rcore

c      Rcore4=Rcore2*Rcore2

c      Rcore16=Rcore4*Rcore4*Rcore4*Rcore4

c      Rcore14=Rcore16/Rcore2

c      cte1=14.d0*Rcore14

c      cte2=1.d0/(rcutoff4*rcutoff4*rcutoff4*rcutoff4)

      polR02=polR0*polR0

      polmaxL2=polmaxL*polmaxL

c      distmax2=(polmaxL+Rcore)**2

      

      cte4=ampli*temp

      cte3=2.d0*cte4/polmaxL2



c      ctep1=(Rcore/rcutoff)**14

c      ctep2=7.d0*ctep1*(Rcore/rcutoff)**2

c      ctep3=(rcutoff/Rcore)**2



c     Now we regularize the polymer-induced force to avoid a numerical instability.

c     This is related to the fact that, due to the discrete character of the time step in

c     the simulation, we can have two particles connected by a polymer at a distance

c     larger than the maximum extension, i.e. L+Rcore. In this case, the force formula

c     we use changes sign, giving rise to an instability. To solve this problem, we

c     regularize the polymer force law by changing it near the divergency. In particular,

c     we introduce an "if" structure, such that whenever (r+Rcore)>fL*L, where fL is 

c     a number smaller but very close to one, e.g. fL=0.999, we substitute the FENE 

c     divergent force by a "ficticious" harmonic force to go back to the physical region. 

c     This force is chosen to be continuous and derivable around the connection point.

c      oort4=1.d0/(fRcore**8)

c      aparam=-cte1*(oort4*oort4-cte2)

      fL2=fL*fL

      fLL=fL*polmaxL

      fLL2=fLL*fLL

      aparam=(-2.d0)*cte4*(1.d0+fL2)/(polmaxL*(1.d0-fL2))**2

      bparam=(-cte3)*polmaxL*fl/(1.d0-fL2) - aparam*fLL

      cparam=cte4*dlog(1.d0/(1.d0-fL2)) +

     &0.5d0*aparam*fLL2 +bparam*fLL



c      write(*,*)aparam,bparam,cparam



c      rsq0=0.99*(polmaxL+Rcore)                                    ! 1% away from the divergence

c      aux0=rsq0-Rcore

c      if(aux0.le.0)then

c         write(*,*)'BAD FICTICIOUS FORCE'

c         stop

c      endif

c      aux02=aux0*aux0/polmaxL2

c      aparam=-cte3*aux0/(rsq0*(1.d0-aux02)*polmaxL)          



c     In order to test the constants, we now write the potential and force functions

c      open(051,file=filepot,status='unknown')

c      dx=0.01d0*(polmaxL+Rcore)

c      rsq0=0.d0

c      do i=1,150

c         rsq0=rsq0+dx

c         rsq=rsq0*rsq0

c         oort3=1.d0/(rsq*rsq*rsq)

c         oort4=1.d0/(rsq*rsq*rsq*rsq)

c         aux=rsq0-Rcore

c         aux2=aux*aux/polmaxL2

c         if( rsq .lt. rcutoff2 ) then                    ! the hard-core interaction is not zero

c            pot=Rcore14*oort3*oort4 - ctep1 + 

c     &ctep2*(rsq/Rcore2 - ctep3 )

c            forc=cte1*(oort4*oort4-cte2)

c         else

c            pot=0

c            forc=0

c         endif

c         if(aux.lt.fLL)then

c            forcepol=(-cte3)*aux/(rsq0*(1.d0-aux2))

c            potpol=(-cte4)*dlog(1.d0-aux2)

c         else

c            forcepol=aparam*aux+bparam

c            potpol= (-0.5d0)*aparam*aux*aux - bparam*aux + cparam

c         endif

c         if(aux.lt.polmaxL)then

c            forcepol=-cte3*aux/(rsq0*(1.d0-aux2))

c            potpol=-cte4*dlog(1.d0-aux2)

c         else

c            forcepol=aparam

c            potpol= - aparam*aux

cc            forcepol=aparam*aux

cc            potpol= - 0.5*aparam*aux*aux

c         endif

cc         potpol=-cte4*dlog(1.d0-aux2)

cc         forcepol=-cte3*aux/(rsq0*(1.d0-aux2))



c         write(051,*)rsq0,pot,potpol,forc,forcepol
c         write(051,*)rsq0,pot,potpol



c      enddo

c      close(051)






c

c we generate the 24 seeds needed by the rng. this is done with the help

c of a bad generator but which is good enough for this purpose. see F.

c James Comp. Phys. Comm. 60, 329 (1990). s(n+1)=mod(23*s(n),10**8+1)

c



c

c to make sure that for different system sizes we get different initial

c configurations we make the seeds dependent on ntot.

c



      iseed=iolseed

      do i=24,1,-1

       iseed=abs(iseed-(npart+npol)*182+iseed/(npart+npol))

       iseed=mod(iseed,100000001)

       iseedn=iseed

       do j=1,22  ! to prevent overflow we add 'by hand'

        iseedn=iseedn+iseed

        if( iseedn .gt. 100000001 ) iseedn=iseedn-100000001

       enddo

       iseed=iseedn

       kseed(i)=abs(mod(iseed,itwop24))

      enddo



c

c we run the rng for a while to get rid of the initial values

c



      do j=1,10000

       iuni=kseed(i24)-kseed(j24)-icarry

       if( iuni .lt. 0 ) then

         iuni=iuni+itwop24

         icarry=1

       else

         icarry=0

       endif

       kseed(i24)=iuni

       i24=i24-1

       if( i24 .eq. 0 ) i24=24

       j24=j24-1

       if( j24 .eq. 0 ) j24=24

      enddo





c

c  read in the old configuration

c



       open(062,file=filein2,status='old',form='formatted')



       enerkin=0.d0

       npp=0

       do i=1,npart


          read(062,*) kk,pos(3*i-2),pos(3*i-1),
     &pos(3*i),vel(3*i-2),vel(3*i-1),vel(3*i)

	write(1009,*)pos(3*i-2),pos(3*i-1),pos(3*i)

          read(062,*) numpol(i)             ! Number of polymers linking particle i to particles j>i (every polymer is counted ONLY ONCE)

          npp=npp+numpol(i)

          enerkin=enerkin+vel(3*i-2)*vel(3*i-2)+vel(3*i-1)*vel(3*i-1)+
     &vel(3*i)*vel(3*i)


          if(numpol(i).ge.idimen)then

             write(*,*)'numpol(i) >=',idimen,' for i=',i

             write(*,*)'PROBLEM WITH VECTORS DIMENSION'

             stop

          endif

          do j=1,numpol(i)                              

             read(062,*)kpol(i,j),neighpol(i,j)           ! kpol(i,j): global index associated to the j-th polymer linked to particle i

             kk=2*kpol(i,j)                    ! neighpol(i,j): Particle connected to particle i by j-th polymer (of a total of numpol(i) polymers)

             kendpoint(kk-1)=i                                    ! kendpoint(kk-1): lower-index particle attached to polymer kpol(i,j)       

             kendpoint(kk)=neighpol(i,j)                      ! kendpoint(kk): higher-index particle attached to polymer kpol(i,j) 

             indexpol(kpol(i,j))=j                                  ! indexpol(kpol(i,j)): local index associated to polymer with global index kpol(i,j) 

          enddo

       enddo

       close(062)



       if(npp.ne.npol)then

          write(*,*)'npp=',npp,'  npol=',npol,

     &'  DIFFERENT!!: CHECK INITIAL CONFIGURATION'

          stop

       endif


c       write(*,*)'enerkin=',0.5d0*enerkin/npart

c       read(*,*)



c

c test whether all particles are in the box

c



      do i=1,npart3

         if( pos(i) .lt. 0 .or. pos(i) .gt. xlbox ) then

            print*,'error1',i,pos(i)

            write(071,*)'#  error1  ',i,pos(i)

            stop

         endif

      enddo

 

      istarttime=0
       

c

c initialize runpos.

c

      do i=1,npart3

         runpos(i)=pos(i)

      enddo



      write(071,*)'#  ',filein2

      call flushkob(071)



c write the first configuration to the dump files



      itime=0

      dump_file=fileout2

      call addnumtostring(dump_file,itime)



      open(075,file=dump_file,status='unknown',form='formatted')

      write(075,*) itime

      do j=1,npart

c         write(075,*) sngl(pos(3*j-2)),sngl(pos(3*j-1)),sngl(pos(3*j))

         write(075,301) j,sngl(pos(3*j-2)),sngl(pos(3*j-1)),

     &sngl(pos(3*j)),sngl(vel(3*j-2)),sngl(vel(3*j-1)),sngl(vel(3*j))

          write(075,*) numpol(j)

          do nothing=1,numpol(j)

             write(075,400) kpol(j,nothing),neighpol(j,nothing)

          enddo         

      enddo



      close(075)







c

c do the md

c



      displamax=0

      acmax=0



      call neilist2(pos)


      call force(pos,acc)

      call analyse1(pos,vel,acc)


      xmed=xmed+1 !what is this ?

      xinsttemp=ekin*2.d0/3.d0

      press=xinsttemp*density+virial*16.d0/xlbox**3

       write(071,fmt='(g12.6,21g23.15)') 

     &  istarttime*hstep,xinsttemp,ekin+epot+epotpol,ekin,epot,epotpol,

     &  press,displamax,acmax,fmax1,fmax2,distm1,distm2,fracloop,0.,0.,

     &  pos(3),pos(6),vel(3),vel(4),vel(5)

      call flushkob(071)



c     Do we have to compute the correlation functions? 

      itime=0


c       if( mod(itime,itimewrite(nwrite)) .eq. 0 ) then
       if( itime.eq.itimewrite(nwrite)) then

c          write(*,*)nwrite,itimewrite(nwrite)

          fichier=filein0

          call addnumtostring(fichier,itime)

          open(075,file=fichier,status='unknown')

          do j=1,npart

             write(075,301) j,sngl(runpos(3*j-2)),sngl(runpos(3*j-1)),

     &sngl(runpos(3*j)),sngl(vel(3*j-2)),sngl(vel(3*j-1)),sngl(vel(3*j))

            write(075,*) numpol(j)

             do nothing=1,numpol(j)

                write(075,400) kpol(j,nothing),neighpol(j,nothing)

             enddo

          enddo

          close(075)

 8989     nwrite=nwrite+1

          if(itimewrite(nwrite).eq.itimewrite(nwrite-1))goto 8989

       endif








c       pause


      do itime=istarttime+1,ntotstep

c         call maxbol(temp,npart,hstep)


c

c update positions and make first step for velocities.

c take into account boundary conditions.

c



c       write(*,*)'itime=',itime,'  xlbox=',xlbox



       do i=1,npart3

          displa=hstep*vel(i)+hsqd2*acc(i)

          pos(i)=pos(i)+displa

          runpos(i)=runpos(i)+displa

          vel(i)=vel(i)+hstepd2*acc(i)

          if( pos(i) .lt. 0    ) then

             pos(i)=pos(i)+xlbox

          else

             if( pos(i) .gt. xlbox ) pos(i)=pos(i)-xlbox

          endif

          displamax=max(abs(displa),displamax)

          acmax=max(abs(acc(i)),acmax)

          

          if(abs(displa).gt.1)write(*,*)'time=',itime*hstep,'  i=',i,

     &'  displa=',displa



       enddo



c       write(34,fmt='(11(1x,f8.3))')itime*hstep,pos(19),pos(20),

c     &pos(21),vel(19),vel(20),vel(21),acc(19),acc(20),acc(21),displamax



c       write(*,*)'time=',itime*hstep,'  displamax=',displamax



       call force(pos,acc)



c

c make second step for the velocities

c



c       ekin=0.d0

       do i=1,npart3

          vel(i)=vel(i)+hstepd2*acc(i)

c          ekin=ekin+vel(i)*vel(i)

       enddo

c       ekin=0.5d0*ekin/npart



c

c now this md step is done. make some analysis or influence the system if needed.

c


c     Do we have to compute the correlation functions? 


c       if(itimewrite(nwrite).gt.0.and.
c     & mod(itime,itimewrite(nwrite)).eq.0)then

	if(itime.eq.itimewrite(nwrite))then

c          write(*,*)nwrite,itimewrite(nwrite),itime


          fichier=filein0

          call addnumtostring(fichier,itime)

          open(075,file=fichier,status='unknown')


          do j=1,npart


             write(075,301) j,sngl(runpos(3*j-2)),sngl(runpos(3*j-1)),

     &sngl(runpos(3*j)),sngl(vel(3*j-2)),sngl(vel(3*j-1)),sngl(vel(3*j))

             write(075,*) numpol(j)

             do nothing=1,numpol(j)

                write(075,400) kpol(j,nothing),neighpol(j,nothing)

             enddo

          enddo

          close(075)

7878     nwrite=nwrite+1

          if(itimewrite(nwrite).eq.itimewrite(nwrite-1))goto 7878

       endif





c do we have to update the neighbor list?



       if( mod(itime,itneighli) .eq. 0 ) call neilist2(pos)
c         call neilist2(pos)





c do we have to thermalize the velocities?



c       if( itime .le. icooltim .and. mod(itime,icoolstep) .eq. 0 ) then


c         call maxbol(temp,npart,hstep)



c         if( itime+icoolstep .gt. icooltim ) then

c this is the last time we cool. let's adjust the total energy if we want to do it

c           if( etot .gt. -10 ) then

c             call analyse1(pos,vel,acc)

c             xmed=xmed+1

c             scafac=sqrt((etot-epot-epotpol)/ekin)

c             do i=1,npart3

c              vel(i)=vel(i)*scafac

c             enddo

c           endif

c         endif

c       endif





c do we have to measure and print the status of the system?



       if(mod(itime,nmedida) .eq. 0) then



         call analyse1(pos,vel,acc)

         xmed=xmed+1



         xinsttemp=ekin*2.d0/3.d0

         press=xinsttemp*density+virial*16.d0/xlbox**3

         faccept=xaccept/xmov



         if(mod(itime,nprintstep) .eq. 0) then



            write(071,fmt='(g12.6,21g23.15)')

     &   itime*hstep,xinsttemp,ekin+epot+epotpol,ekin,epot,epotpol,

     &   press,displamax,acmax,fmax1,fmax2,distm1,distm2,fracloop,

     &   faccept,xcte/faccept,pos(3),pos(6),vel(3),vel(4),vel(5)

            call flushkob(071)



         endif



       endif





c do we have to write the configuration?


c       if(mod(itime,ndumpstep) .eq. 0 .and. fileout1 .ne. ' ') then

       if(mod(itime,ndumpstep) .eq. 0 ) then

        

         dump_file=fileout2

         call addnumtostring(dump_file,itime)



         open(075,file=dump_file,status='unknown',form='formatted')

         write(075,*) itime

         do j=1,npart

            write(075,301) j,sngl(pos(3*j-2)),sngl(pos(3*j-1)),

     &sngl(pos(3*j)),sngl(vel(3*j-2)),sngl(vel(3*j-1)),sngl(vel(3*j))

            write(075,*) numpol(j)

            do nothing=1,numpol(j)

               write(075,400) kpol(j,nothing),neighpol(j,nothing)

            enddo

         enddo

         close(075)







       endif



c

c test whether all particles are still in the box

c



       do i=1,npart3

        if( pos(i) .lt. 0 .or. pos(i) .gt. xlbox ) then

          print*,'error2',i,pos(i)

          write(071,*)'#  error2   ',i,pos(i)

          stop

        endif

       enddo

c      write(800,*)itime, virial

      enddo                    ! time loop





c

c write the file with the final state

c



      if( fileout3 .ne. ' ' ) then

         open(072,file=fileout3,status='unknown',form='formatted')

c        do i=1,ntot

c         write(072,fmt='(6g25.15)') pos(3*i-2),pos(3*i-1),pos(3*i),

c     &                vel(3*i-2),vel(3*i-1),vel(3*i)

c        enddo

         do i=1,npart

            write(072,301) i,sngl(pos(3*i-2)),sngl(pos(3*i-1)),

     &sngl(pos(3*i)),sngl(vel(3*i-2)),sngl(vel(3*i-1)),sngl(vel(3*i))

            write(072,*) numpol(i)

            do nothing=1,numpol(i)

               write(072,400) kpol(i,nothing),neighpol(i,nothing)

            enddo

         enddo

         close(072)

      endif





      print*,'good exit'

      write(071,*)'#  good exit'

      close(071)

  

      close(088)



      stop



 100  format(5(1x,e25.17))

 101  format(3(1x,e25.17))

 301  format(1x,i6.6,6(1x,e14.7))

 400  format(2(1x,i6.6))



      end





 

      subroutine force(pos,acc)



      implicit double precision (a-h,o-z)



      parameter (ipartdim1=10000,ipartdim=3*ipartdim1,idimen=600,

     &idimen2=150)



      dimension neighbor(idimen2*ipartdim1),neighpol(ipartdim1,idimen),

     &numpol(ipartdim1)

      dimension pos(ipartdim),acc(ipartdim)
 
      dimension dia(ipartdim1)


      common/kob1/ xlbox,xlboxd2,npart,npol,ntot

      common/kob4/ rcutoff, rcutoff2

      common/kob6/ neighbor

      common/kob11/ virial

      common/pablo1/ npart3

      common/pablo2/ neighpol,numpol

      common/pablo3/ cte1,cte2,cte3,cte4,polmaxL,polmaxL2,Rcore,epsilon

      common/pablo7/ aparam,bparam,cparam,fLL,fLL2

      common/pablo8/ hkk,fmax1,fmax2,distm1,distm2,itime

      common/pc1/dia



      virial=0.d0


      rcutoff4=rcutoff2*rcutoff2


c

c first zero all accelerations

c



      do i=1,npart3

         acc(i)=0.d0

      enddo

c Now get the fluid velocity at the position of the particles 
c      call get_fluid_velocity(pos,fluid_vel)

c

c compute the forces due to the hard-core and the polymer-induced interactions between particles

c



      ihigh=0

      do i=1,npart

         xx=pos(3*i-2)

         yy=pos(3*i-1)

         zz=pos(3*i)

         fxi=0.d0

         fyi=0.d0

         fzi=0.d0

c

c     First, hard-core interactions

c

         low=ihigh+2                                 ! Position of the first neighbor of particle i in vector neighbor().

         ihigh=ihigh+1+neighbor(low-1)    ! Position of the last neighbor of particle i in vector neighbor().



         do neighno=low,ihigh

            j=neighbor(neighno)

            dx=xx-pos(3*j-2)

            dy=yy-pos(3*j-1)

            dz=zz-pos(3*j)



            if( dx .gt. xlboxd2 ) then 

               dx=dx-xlbox

            else

               if( dx .lt. -xlboxd2 ) dx=dx+xlbox

            endif



            if( dy .gt. xlboxd2 ) then

               dy=dy-xlbox

            else

               if( dy .lt. -xlboxd2 ) dy=dy+xlbox

            endif



            if( dz .gt. xlboxd2 ) then

               dz=dz-xlbox

            else

               if( dz .lt. -xlboxd2 ) dz=dz+xlbox

            endif

            rsq=0.0d0
            rsq=dx*dx+dy*dy+dz*dz         ! Center-to-center distance between particles

             Rcore=0.0d0
             Rcore=0.50d0*(dia(i)+dia(j))

             Rcore2=Rcore*Rcore

             Rcore4=Rcore2*Rcore2

             Rcore16=Rcore4*Rcore4*Rcore4*Rcore4

             Rcore14=Rcore16/Rcore2

            rcutoffp4=0.0d0
            rcutoffp4=rcutoff4*Rcore4

            cte1=14.d0*Rcore14
            cte2=1.d0/(rcutoffp4*rcutoffp4*rcutoffp4*rcutoffp4)


            rcutoffp2=0.0d0
            rcutoffp2=rcutoff2*Rcore2

            if( rsq .lt. rcutoffp2 ) then             ! the interaction is not zero

               oort4=1.d0/(rsq*rsq*rsq*rsq)

               tforce=cte1*(oort4*oort4-cte2)

               virial=virial+tforce*rsq

               fxi=fxi+dx*tforce                              ! REPULSIVE FORCE: (+) SIGN

               acc(3*j-2)=acc(3*j-2)-dx*tforce     ! Force(x) of i particle on j particle

               

               fyi=fyi+dy*tforce                              ! REPULSIVE FORCE: (+) SIGN

               acc(3*j-1)=acc(3*j-1)-dy*tforce     ! Force(y) of i particle on j particle

               

               fzi=fzi+dz*tforce                              ! REPULSIVE FORCE: (+) SIGN

               acc(3*j)=acc(3*j)-dz*tforce           ! Force(z) of i particle on j particle



               fmax1=max(fmax1,dabs(tforce))

               distm1=min(distm1,dsqrt(rsq))



c               if(i.eq.7)then

c                  write(35,*)'C: time=',itime*hstep,'  i=',i,'  j=',

c     &j,'  r=',dsqrt(rsq),'  Fx=',dx*tforce,'  Fy=',dy*tforce,

c     &'  Fz=',dz*tforce                  

c               elseif(j.eq.7)then

c                  write(35,*)'C: time=',itime*hstep,'  i=',i,'  j=',

c     &j,'  r=',dsqrt(rsq),'  Fx=',-dx*tforce,'  Fy=',-dy*tforce,

c     &'  Fz=',-dz*tforce

c               endif



            endif




         enddo                  ! Loop over hard-core neighbors




         acc(3*i-2)=acc(3*i-2)+fxi                     ! Force(x) of all neighbor particles on i particle

         acc(3*i-1)=acc(3*i-1)+fyi                     ! Force(y) of all neighbor particles on i particle

         acc(3*i)=acc(3*i)+fzi                           ! Force(z) of all neighbor particles on i particle

c Now add the drag force. 

c               acc(3*j-2)= tau_drag1*(vel(3*j-2)  - fluid_vel(3*j-2))
c              acc(3*j-1)= tau_drag1*(vel(3*j-1)  - fluid_vel(3*j-1))
c               acc(3*j-1)= tau_drag1*(vel(3*j-2)  - fluid_vel(3*j-1))



         

      enddo                     ! over particles



c      write(35,*)'%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%'






      return

      end







      subroutine analyse1(pos,vel,acc)



      implicit double precision (a-h,o-z)



      parameter (ipartdim1=10000,ipartdim=3*ipartdim1,idimen=600,

     &idimen2=150)



      dimension neighbor(idimen2*ipartdim1),numpol(ipartdim1),

     &neighpol(ipartdim1,idimen)

      dimension pos(ipartdim),vel(ipartdim),acc(ipartdim)

      dimension  histox(-100:100),histoy(-100:100),histoz(-100:100),

     &histor(0:100),histovx(-100:100),histovy(-100:100),

     &histovz(-100:100),histov(0:100)

      dimension dia(ipartdim1)


      common/kob1/ xlbox,xlboxd2,npart,npol,ntot

      common/kob3/ hstep,ekin,epot,epotpol,fracloop

      common/kob4/ rcutoff, rcutoff2

      common/kob6/ neighbor

      common/pablo1/ npart3

      common/pablo2/ neighpol,numpol

      common/pablo3/ cte1,cte2,cte3,cte4,polmaxL,polmaxL2,Rcore,epsilon

      common/pablo4/ ctep1,ctep2,ctep3,Rcore2,Rcore14

      common/pablo7/ aparam,bparam,cparam,fLL,fLL2

      common/measure1/ histox,histoy,histoz,histor

      common/measure2/ histovx,histovy,histovz,histov

c      common/measure3/ histoE,histoEv,hstoEc,histoEp

      common/mpara/ vmax,dhv,dhx,xpmed,npoints

      common/pc1/dia

c

c compute kinetic energy

c

      ekin=0

      do i=1,npart3

       ekin=ekin+vel(i)*vel(i)

      enddo

      ekin=0.5d0*ekin/npart

c

c now compute the potential energy

c



c      do i=1,npart

c         write(*,*)'i=',i,'  numpol(i)=',numpol(i)

c         do j=1,numpol(i)

c            write(*,*)j,neighpol(i,j)

c         enddo

c      enddo

c      read(*,*)


      epot=0

      epotpol=0

      ihigh=0

      nlooped=0

      do i=1,npart


         xx=pos(3*i-2)

         yy=pos(3*i-1)

         zz=pos(3*i)

         low=ihigh+2

         ihigh=ihigh+1+neighbor(low-1)



         do neighno=low,ihigh

            j=neighbor(neighno)

            dx=xx-pos(3*j-2)

            dy=yy-pos(3*j-1)

            dz=zz-pos(3*j)



            if( dx .lt. -xlboxd2 ) dx=dx+xlbox

            if( dx .gt.  xlboxd2 ) dx=dx-xlbox

            if( dy .lt. -xlboxd2 ) dy=dy+xlbox

            if( dy .gt.  xlboxd2 ) dy=dy-xlbox

            if( dz .lt. -xlboxd2 ) dz=dz+xlbox

            if( dz .gt.  xlboxd2 ) dz=dz-xlbox


            rsq=0.0d0
            rsq=dx*dx+dy*dy+dz*dz

            oort3=1.d0/(rsq*rsq*rsq)

            Rcore=0.0d0
            Rcore=0.50d0*(dia(i)+dia(j))

            Rcore2=Rcore*Rcore
	   
            Rcore4=Rcore2*Rcore2

            Rcore16=Rcore4*Rcore4*Rcore4*Rcore4

            Rcore14=Rcore16/Rcore2

            rcutoffp=rcutoff*Rcore
            rcutoffp2=rcutoff2*Rcore2

               ctep1=(Rcore/rcutoffp)**14

               ctep2=7.d0*ctep1*(Rcore/rcutoffp)**2

               ctep3=(rcutoffp/Rcore)**2
 



            if( rsq .lt. rcutoffp2 ) then ! the interaction is not zero

               epot=epot + Rcore14*oort3*oort3/rsq - ctep1 + 

     &ctep2*(rsq/Rcore2 - ctep3 )



            endif



c            write(35,*)dsqrt(rsq),rsq,epot



         enddo                  ! hard-core neighbors




      enddo                     ! particles



      epot=epot/npart

      epotpol=epotpol/npart

      fracloop=dfloat(nlooped)/npol

c	print*,epot, epotpol

      return

      end









      subroutine neilist2(pos)



      implicit double precision (a-h,o-z)



      parameter (ipartdim1=10000,ipartdim=3*ipartdim1)

      parameter (idimen=600,ipoldim=20*ipartdim1,idimen2=150)



      logical*1 first



      dimension neighbor(idimen2*ipartdim1),

     &neighboro(idimen2*ipartdim1),nvecp(ipartdim1),

     &kvecpolymer(ipartdim1,idimen)

      dimension pos(ipartdim)

      dimension dia(ipartdim1)



      common/kob1/ xlbox,xlboxd2,npart,npol,ntot

      common/kob4/ rcutoff, rcutoff2

      common/kob6/ neighbor

      common/kob7/ rskin2

      common/kob10/ neighboro

      common/pablo5/ distmax2,kcountmax

      common/pablo6/ nvecp,kvecpolymer
    
      common/pablo3/ cte1,cte2,cte3,cte4,polmaxL,polmaxL2,Rcore,epsilon
     

	common/pc1/dia



      data first /.true./





      do i=1,npart

         nvecp(i)=0               ! Number of polymer-mediated neighboring particles of particle i  

      enddo



      noneigh=1

      kcount=1

      ihigh=0

      do i=1,npart

       xx=pos(3*i-2)

       yy=pos(3*i-1)

       zz=pos(3*i)

       low=ihigh+2

       ihigh=ihigh+1+neighboro(low-1)



c     Since we can create polymer loops, linking one particle with itself, we have to count 

c     particle i as one of its own polymer-mediated neighbor. In this way, the first polymer-mediated

c     neighbor particle is always the given particle.

       nvecp(i)=nvecp(i)+1

       kvecpolymer(i,nvecp(i))=i



       do j=i+1,npart               ! Every pair of interacting particles is counted ONLY ONCE.

           
c          Rcore=1.0d0
c --------------------------------------------------------
           Rcore=0.0d0
           Rcore=0.50d0*(dia(i)+dia(j))

          distmax2=0.0d0
          distmax2=(polmaxL+Rcore)**2

          rcutoffp2=0.0d0
          rcutoffp2=rcutoff2*Rcore*Rcore

           rskinp2=0.0d0
           rskinp2=rskin2*Rcore*Rcore
c ---------------------------------------------------------

          dx=xx-pos(3*j-2)


          if( dx .lt. -xlboxd2 ) then

             dx=dx+xlbox

          else

             if( dx .gt.  xlboxd2 ) dx=dx-xlbox

          endif



c        if( dx*dx .ge. rskinp2 ) goto 12

          if( dx*dx .ge. distmax2 ) goto 122



          dy=yy-pos(3*j-1)

          if( dy .lt. -xlboxd2 ) then

             dy=dy+xlbox

          else

             if( dy .gt.  xlboxd2 ) dy=dy-xlbox

          endif



          rsq=dx*dx+dy*dy

c        if( rsq .ge. rskinp2 ) goto 12

          if( rsq .ge. distmax2 ) goto 122



          dz=zz-pos(3*j)

          if( dz .lt. -xlboxd2 ) then

             dz=dz+xlbox

          else

             if( dz .gt.  xlboxd2 ) dz=dz-xlbox

          endif


          rsq=rsq+dz*dz


          if( rsq .le. distmax2 ) then                  ! we found a polymer-mediated neighbor 



             nvecp(i)=nvecp(i)+1

             nvecp(j)=nvecp(j)+1



             kvecpolymer(i,nvecp(i))=j

             kvecpolymer(j,nvecp(j))=i

           

             if( rsq .lt. rskinp2 ) then                     ! we found a Verlet neighbor


                kcount=kcount+1 

                neighbor(kcount)=j

c

c test whether the particle sneaked into the interaction range

c

          if( first ) goto 312


                if( rsq .lt. rcutoffp2 ) then

                   do k=low,ihigh

                      if( neighboro(k) .eq. j ) goto 312

                   enddo

                   print*,'error4',i,j,'  first=',first,'  rsq0=',

     &dsqrt(rsq)

                   write(071,*)'#  error4  ',i,j,'  first=',first

                endif

 312         endif



          endif



 122   enddo



       neighbor(noneigh)=kcount-noneigh

       kcount=kcount+1

       noneigh=kcount

       

      enddo



c      if( kcount .gt. 150000 ) then

      if( kcount .gt. kcountmax ) then

        print*,'error3', kcount

        write(071,*)'#  error3  ',kcount

      endif



      do i=1,kcount

       neighboro(i)=neighbor(i)

      enddo



      first=.false.


      return

      end





c

c generate a maxwell-boltzmann distribution for the velocities: WITH ZERO TOTAL MOMENTUM 

c

      subroutine maxbol(temp,npart,hstep)



      implicit double precision (a-h,o-z)



      parameter (ipartdim1=10000,ipartdim=3*ipartdim1)



      dimension kseed(24)

      dimension vel(ipartdim)



      common/kob8/ vel

      common/kob9/ tdtp24,icarry,i24,j24,kseed,itwop24

      common/pablo1/ npart3



      do i=1,npart3,2

   10  iuni=kseed(i24)-kseed(j24)-icarry

       if( iuni .lt. 0 ) then

         iuni=iuni+itwop24

         icarry=1

       else

         icarry=0

       endif

       kseed(i24)=iuni

       i24=i24-1

       if( i24 .eq. 0 ) i24=24

       j24=j24-1

       if( j24 .eq. 0 ) j24=24

       v1=tdtp24*(iuni+1)-1.0



       iuni=kseed(i24)-kseed(j24)-icarry

       if( iuni .lt. 0 ) then

         iuni=iuni+itwop24

         icarry=1

       else

         icarry=0

        endif

       kseed(i24)=iuni

       i24=i24-1

       if( i24 .eq. 0 ) i24=24

       j24=j24-1

       if( j24 .eq. 0 ) j24=24

       v2=tdtp24*(iuni+1)-1.0



       s=v1*v1+v2*v2

       if( s .ge. 1.0 ) goto 10

       r=sqrt((-2.0)*log(s)/s)



       vel(i)=v1*r

       vel(i+1)=v2*r



      enddo



      cmvx=0.d0

      cmvy=0.d0

      cmvz=0.d0

      do i=1,npart

       cmvx=cmvx+vel(3*i-2)

       cmvy=cmvy+vel(3*i-1)

       cmvz=cmvz+vel(3*i)

      enddo



      cmvx=cmvx/npart

      cmvy=cmvy/npart

      cmvz=cmvz/npart



c     Ensure zero total momentum.

      do i=1,npart

       vel(3*i-2)=vel(3*i-2)-cmvx

       vel(3*i-1)=vel(3*i-1)-cmvy

       vel(3*i  )=vel(3*i  )-cmvz

      enddo



      ekin=0.d0

      do i=1,npart3

       ekin=ekin+vel(i)*vel(i)

      enddo



c      write(*,*)'Energy antes=',0.5d0*ekin/ntot



c      scafac=sqrt(temp*ntot/ekin/16.d0)

      scafac=sqrt(3.d0*temp*npart/ekin)

c	print*,ekin,scafac



      ekin2=0.d0

      do i=1,npart3

         vel(i)=vel(i)*scafac

c       print*,vel(i)

         ekin2=ekin2+vel(i)*vel(i)

      enddo



c      write(*,*)'Energy despues=',0.5d0*ekin2/npart,

c     &'   E equipart=',1.5d0*temp

c      read(*,*)

c	print*,ekin/npart


      return

      end



cc

cc

cc this subroutine forces the contents of the iobuffer to be

cc writen to the file

cc

cc





      subroutine flushkob(iounit)



      implicit double precision (a-h,o-z)



      character*1 chara



      endfile iounit

      rewind iounit

10    read (iounit,fmt='(a)', end=40) chara

      goto 10

40    return

      end





cc

cc

cc this subroutine adds the number 'number' to the string 'string'

cc

cc



      subroutine addnumtostring(string,number)



      implicit double precision (a-h,o-z)



      character*(*) string



      isnum=number

      do i=len(string),1,-1

       if (string(i:i) .ne. ' ') goto 10

      enddo

   10 istrlen=i





      nodig=int(log10(1.0*isnum+0.2))+1

      do i=nodig,1,-1

       num=isnum/10**(i-1)

       string(istrlen+1:istrlen+1)=char(48+num)

       istrlen=istrlen+1

       isnum=isnum-num*10**(i-1)

      enddo



      return



      end





      SUBROUTINE sort(n,arr)

      INTEGER n,M,NSTACK

c      REAL arr(n)

      INTEGER arr(n),a,temp

      PARAMETER (M=7,NSTACK=50)

      INTEGER i,ir,j,jstack,k,l,istack(NSTACK)

c      REAL a,temp

      jstack=0

      l=1

      ir=n

1     if(ir-l.lt.M)then

        do 12 j=l+1,ir

          a=arr(j)

          do 11 i=j-1,1,-1

            if(arr(i).le.a)goto 2

            arr(i+1)=arr(i)

11        continue

          i=0

2         arr(i+1)=a

12      continue

        if(jstack.eq.0)return

        ir=istack(jstack)

        l=istack(jstack-1)

        jstack=jstack-2

      else

        k=(l+ir)/2

        temp=arr(k)

        arr(k)=arr(l+1)

        arr(l+1)=temp

        if(arr(l+1).gt.arr(ir))then

          temp=arr(l+1)

          arr(l+1)=arr(ir)

          arr(ir)=temp

        endif

        if(arr(l).gt.arr(ir))then

          temp=arr(l)

          arr(l)=arr(ir)

          arr(ir)=temp

        endif

        if(arr(l+1).gt.arr(l))then

          temp=arr(l+1)

          arr(l+1)=arr(l)

          arr(l)=temp

        endif

        i=l+1

        j=ir

        a=arr(l)

3       continue

          i=i+1

        if(arr(i).lt.a)goto 3

4       continue

          j=j-1

        if(arr(j).gt.a)goto 4

        if(j.lt.i)goto 5

        temp=arr(i)

        arr(i)=arr(j)

        arr(j)=temp

        goto 3

5       arr(l)=arr(j)

        arr(j)=a

        jstack=jstack+2

        if(jstack.gt.NSTACK)pause 'NSTACK too small in sort'

        if(ir-i+1.ge.j-l)then

          istack(jstack)=ir

          istack(jstack-1)=i

          ir=j-1

        else

          istack(jstack)=j-1

          istack(jstack-1)=l

          l=i

        endif

      endif

      goto 1

      END

